
MODULE module_chemistry_driver

    USE rconstants , ONLY: &
         pi180            ! (IN)

    USE mem_grid, ONLY: &
         time,          & ! (IN)
         iyear1,        & ! (IN)
         imonth1,       & ! (IN)
         idate1,        & ! (IN)
         itime1,        & ! (IN)
         dtlt,          & ! (IN)
         ngrid,         & ! (IN)
         ngrids,        & ! (IN)
         dtlongn,       & ! (IN)
         dzt,           & ! (IN)
         zm,            & ! (IN)
         initial,       & ! (IN)
         zt,            & ! (IN)
         centlat,       & ! (IN)
         centlon,       & ! (IN)
         grid_g,        & ! %glat(IN), %glon(IN), %rtgt(IN), %lpw(IN), %dxt(IN)
         platn,         & ! (IN)
         plonn,         & ! (IN)
         deltaxn,       & ! (IN)
         deltayn,       & ! (IN)
         xt,            & ! (IN)
         yt               ! (IN)

    USE chem1_list, ONLY :          &
         chem_nspecies =>nspecies,  & ! (IN)
         spc_chem_alloc=>spc_alloc, & ! (IN)
         spc_chem_name =>spc_name,  & ! (IN)
         src,                       & ! (IN)
         on,                        & ! (IN)
         off,                       & ! (IN)
         transport,                 & ! (IN)
         nspecies,                  & ! (IN)
         spc_alloc,                 & ! (IN)
         no,                        & ! (IN)
         weight,                    & ! (IN)
         co,                        & ! (IN)
         photojmethod,              & ! (IN)
         nr,                        & ! (IN)
         nr_photo                     ! (IN)

    USE mem_chem1, ONLY:               &
         CHEMISTRY,                    & ! (IN)
         maxnspecies,                  & ! (IN)
         nspecies_chem_transported,    & ! (IN)
         nspecies_chem_no_transported, & ! (IN)
         transp_chem_index,            & ! (IN)
         no_transp_chem_index,         & ! (IN)
         nsrc,                         & ! (IN)
         src_name,                     & ! (IN)
         antro,                        & ! (IN)
         bburn,                        & ! (IN)
         bioge,                        & ! (IN)
         chem1_src_z_dim_g,            & ! (IN)
         chem1_g,                      & ! %sc_t(INOUT), %sc_p(INOUT)
         chem1_src_g,                  & ! %sc_src(INOUT)
         chem1_src_vars,               & ! Type
         chem1_vars,                   & ! Type
         N_DYN_CHEM_N,                 & ! (IN)
         N_DYN_CHEM,                   & ! (INOUT)
         split_method,                 & ! (IN)
         isplit                          ! (IN)


    USE mem_chem1aq, ONLY: &
         CHEMISTRY_AQ        ! (IN)

    USE mem_basic, ONLY: &
         basic_g,        &   ! %up(IN), %vp(IN), %wp(IN), %theta(IN), %pp(IN), %pi0(IN), %dn0(IN?), %rv(IN), %rtp(IN)
         basic_vars          ! Type

!    USE chem_sources, ONLY: &
!         plumerise,         & ! (IN)
!         pr_time=>PRFRQ       ! (IN)

    USE mem_radiate, ONLY : &
         radfrq,            & ! (IN)
         lonrad,            & ! (IN) (nao usado)
         radiate_g,         & ! %cosz(IN), %rlongup(IN), %albedt(IN)
         radiate_vars,      & ! Type
         ilwrtyp,           & ! (IN)
         iswrtyp              ! (IN)

    USE FastJX, ONLY :  &
         FastJX_driver, & ! Subroutine
         fast_JX_g        ! %jphoto(IN)

    USE node_mod, ONLY : &
         mynum,          & ! (IN)
         mmxp,           & ! (IN)
         mmyp,           & ! (IN)
         mmzp,           & ! (IN)
         mia,            & ! (IN)
         miz,            & ! (IN)
         mja,            & ! (IN)
         mjz,            & ! (IN)
         mibcon,         & ! (IN)
         mi0,            & ! (IN)
         mj0,            & ! (IN)
         maxgrds,        & ! (IN)
         i0,             & ! (IN)
         j0,             & ! (IN)
         ibcon             ! (IN)

    USE mem_micro, ONLY: &
         micro_g,        & ! %rcp(IN)
         micro_vars        ! Type

    USE micphys, ONLY: &
         level           ! (IN)

    USE grid_dims, ONLY: &
         nzpmax,         &  ! (IN)
         nxpmax,         &  ! (IN)
         nypmax             ! (IN)

    USE mem_scratch1_grell, ONLY: &
         ierr4d                     ! (IN)

    USE mem_grell_param, ONLY: &
         mgmxp,                & ! (IN)
         mgmyp,                & ! (IN)
         maxiens,              & ! (IN)
         ngrids_cp               ! (IN)

    USE mem_cuparm, ONLY: &
         nnqparm            ! (IN)

    USE mem_stilt, ONLY: &
         iexev,          & ! (IN)
         stilt_g           ! %dnp(IN)

    USE mem_globrad, ONLY: &
         raddatfn            ! (IN)

    USE carma_fastjx, ONLY: &
         do3,               & ! (IN)
         daer,              & ! (IN)
         na                   ! (IN)

    USE aer1_list, ONLY :          &
         aer_nspecies=>nspecies,   & ! (IN)
         spc_aer_alloc=>spc_alloc, & ! (IN)
         nmodes,                   & ! (IN)
         spc_aer_name =>aer_name,  & ! (IN)
         aer_bburn => bburn,       & ! (IN)
         aer_sdust => sdust,       & ! (IN)
         aer_urban => urban,       & ! (IN)
         aer_bioge => bioge,       & ! (IN)
         aer_marin => marin          ! (IN)

    USE mem_aer1, ONLY:                   &
         aer_nvert_src=>aer1_src_z_dim_g, & ! (IN)
         AEROSOL,                         & ! (IN)
         aer1_g,                          & ! %sc_src(INOUT)
         aer1_vars                          ! Type

!    USE mem_plume_chem1, ONLY: &
!         nveg_agreg,           & ! (IN)
!         tropical_forest,      & ! (IN)
!         boreal_forest,        & ! (IN)
!         savannah,             & ! (IN)
!         grassland,            & ! (IN)
!         plume_mean_g,         & ! %flam_frac(IN), %fire_size(IN)
!         plume_mean_vars         ! Type


!    USE plumegen_coms,ONLY:   &
!         nkp,                 & ! (IN)
!         ntime,               & ! (IN)
!         zt_plumegen=>zt,     & ! (INOUT)
!         zm_plumegen=>zm,     & ! (INOUT)
!         dzt_plumegen=>dzt,   & ! (INOUT)
!         dzm,                 & ! (INOUT)
!         dz,                  & ! (INOUT)
!         ncall,               & ! (INOUT)
!         nrec,                & ! (INOUT)
!         n_setgrid              ! (INOUT)

    USE spack_utils, ONLY: &
         allocindex,       & ! Subroutine
         index_g,          & ! IN
         nob,              & ! IN
         maxblock_size       ! IN

    USE mem_spack, ONLY: &
         alloc_spack,    & ! Subroutine
         spack_alloc,    & ! IN
         spack,          & ! OUT
         spack_2d,       & ! INOUT
         atol,           & ! IN
         rtol              ! IN

    USE solve_sparse, ONLY: &
         get_non_zeros

    USE mod_chem_spack_ros_dyndt, ONLY: &
         chem_ros_dyndt    ! Subroutine

    USE mod_chem_spack_ros, ONLY: &
         chem_ros          ! Subroutine

    USE mod_chem_spack_qssa, ONLY: &
         chem_qssa         ! Subroutine

    USE mod_chem_trans_gasaq, ONLY: &
         trans_gaq         ! Subroutine

    USE mod_chem_orage, ONLY: &
         eclair_driver     ! Subroutine

    USE mod_chem_spack_rodas3_dyndt, ONLY: &
         chem_rodas3_dyndt ! Subroutine

    USE uv_atten, ONLY: &
         uv_attenuation    ! Subroutine

    USE mod_chem_trans_liq, ONLY: &
         trans_liq         ! Subroutine

    USE chem1aq_list, ONLY: &
         nspeciesaq,        & ! (IN)
         ind_gas              ! (IN)

    USE mem_chemic, ONLY: &
         chemic_g            ! (IN)

    USE mem_chem1aq, ONLY: &
         chem1aq_g           ! (INOUT)

    USE mem_aerad, ONLY: &
         nwave              ! (IN)

    USE mem_carma, ONLY: &
         carma              ! (IN)
       
!lfr - tuv
    USE ModtuvDriver, ONLY: tuvDriver
    USE ModTuv      , ONLY: kz,kw,kt,ks,kj
    use tuvParameter, ONLY: nstr,wstart,wstop,nwint,listFiles
!lfr - tuv

    IMPLICIT NONE

    PRIVATE

    PUBLIC ::  chemistry_driver ! Subroutine


CONTAINS

  !------------------------------------------------------------------
  SUBROUTINE chemistry_driver(mzp,mxp,myp,ia,iz,ja,jz,task,nt)

    INTEGER , INTENT(IN) :: mzp
    INTEGER , INTENT(IN) :: mxp
    INTEGER , INTENT(IN) :: myp
    INTEGER , INTENT(IN) :: ia
    INTEGER , INTENT(IN) :: iz
    INTEGER , INTENT(IN) :: ja
    INTEGER , INTENT(IN) :: jz
    INTEGER , INTENT(IN) :: task
    INTEGER , INTENT(IN) :: nt

    INTEGER, DIMENSION(ngrids) :: itchim
    INTEGER, DIMENSION(ngrids) :: ncycle

    INTEGER,PARAMETER :: block_size_qssa=45
    INTEGER,PARAMETER :: block_size_ros=2*32

!lfr - tuv
     !***********************************************************************
     !For test tuv_driver:  
     INTEGER :: nnj
     CHARACTER(len=50) :: sslabel(ks)  ! Spectral weighting functions labels
     CHARACTER(len=50) :: jjlabel(kj)  ! Photolysis coefficients labels
     REAL :: rate(ia:iz,ja:jz,ks,kz)   ! Weighted irradiances (dose rates) W m-2
     REAL :: valj(ia:iz,ja:jz,kj,kz)   ! Photolysis coefficients (j-values)
     REAL :: sirrad(ia:iz,ja:jz,kz,kw) ! Spectral irradiance, [W m-2 nm-1]
     REAL :: saflux(ia:iz,ja:jz,kz,kw) ! Spectral actinic flux, quanta s-1 nm-1 cm-2
     INTEGER,PARAMETER :: blockSize=1
     INTEGER :: i,j,k,jn
  !***********************************************************************				
!lfr - tuv
   !print*,'81'; call flush(6)	  


    !-  call NO production by "eclair" 
    IF( task == 2 .AND. CHEMISTRY >= 0) THEN

       IF(ASSOCIATED(chem1_g(NO,ngrid)%sc_t)) THEN

          CALL eclair_driver(mzp,mxp,myp,ia,iz,ja,jz,maxgrds,dtlt,ngrid,grid_g(ngrid)%rtgt, &
                             grid_g(ngrid)%dxt,grid_g(ngrid)%dyt,dzt,zt,nzpmax,mgmxp,mgmyp, &
                             maxiens,ngrids_cp,ierr4d,nnqparm,basic_g(ngrid)%wp, &
                             basic_g(ngrid)%rtp,basic_g(ngrid)%rv,basic_g(ngrid)%pp, &
                             basic_g(ngrid)%pi0,basic_g(ngrid)%theta,chem1_g(NO,ngrid)%sc_t, &
                             weight,nspecies,no)

       ENDIF

       RETURN

    ELSEIF( task == 3 .AND. CHEMISTRY >= 0) THEN


       IF (CHEMISTRY == 0) THEN
             
          CALL include_tracer_lifetime(mzp,mxp,myp,spc_alloc,nspecies,transport,on,co, &
                                       chem1_g(CO,ngrid)%sc_p(:,:,:), chem1_g(CO,ngrid)%sc_t(:))
      	   
       ELSE

          !-  call chemistry production/loss
           	  
          !- Determining number of blocks and masks, allocating spack type
          IF(.NOT. spack_alloc) THEN 
             CALL allocIndex(block_size_ros,mmxp,mmyp,mmzp,mia,miz,mja,mjz,mibcon,&
                             ngrids,dtlongn,N_DYN_CHEM) 
             CALL alloc_spack(CHEMISTRY)
          END IF
		   
          !- photolysis or attenuation calculations			 
          IF (MOD(time + .001,radfrq) .LT. dtlt .OR. time .LT. 0.001) THEN 
             IF (TRIM(PhotojMethod) == 'FAST-JX') THEN

                CALL FastJX_driver(mzp,mxp,myp,ia,iz,ja,jz,nzpmax,ngrid,ngrids,grid_g(ngrid)%lpw, &
                                   grid_g(ngrid)%rtgt,zm,imonth1,idate1,iyear1,time,dzt,mmzp,mmxp, &
           	  	  	   mmyp,i0,j0,maxgrds,basic_g(ngrid)%pp,basic_g(ngrid)%pi0, &
           		           basic_g(ngrid)%theta,basic_g(ngrid)%dn0,radiate_g(ngrid)%rlongup, &
           		  	   radiate_g(ngrid)%cosz,radiate_g(ngrid)%albedt,raddatfn,do3,daer,na)
	   
             ELSEIF(TRIM(PhotojMethod) == 'LUT') THEN
  !print*,'82'; call flush(6)	  
	       
                CALL uv_attenuation(ngrid,mxp,myp,ia,iz,ja,jz,i0,j0,ngrids,platn,plonn,deltaxn, &
                                    deltayn,nxpmax,nypmax,xt,yt,mmxp,mmyp,nwave,carma(ngrid)%aot, &
                                    ilwrtyp,iswrtyp)
	    
  !print*,'83'; call flush(6)	  
	    ELSEIF(TRIM(PhotojMethod) == 'FAST-TUV') THEN
                
	       call tuvDriver(mzp,mxp,myp,ia,iz,ja,jz,nstr,wstart,wstop,nwint,nnj,jjlabel,sslabel,&
	                      rate,valj,sirrad,saflux,blockSize,listFiles)
	     
             ELSE
	        STOP ' unknonw photolysis calculation method '
           	      
             ENDIF
          ENDIF
             
          !-srf: including options for splliting operator
          !- current dynamic splliting
          N_DYN_CHEM=N_DYN_CHEM_N(ngrid)
	   	     
          IF(split_method == 'PARALLEL' .AND. N_DYN_CHEM > 1) &
               CALL chem_accum(ngrid,mzp,mxp,myp,ia,iz,ja,jz,1,nspecies_chem_transported,chem1_g, &
                               transp_chem_index,N_DYN_CHEM)
          
          IF (MOD(time-(N_DYN_CHEM/ISPLIT-1)*dtlt , N_DYN_CHEM*dtlt) == 0. ) THEN
	   
 	   
             IF(CHEMISTRY == 1) THEN  
                itchim(ngrid)= 75  ! teste 60
                !STOP 'for now use only CHEMISTRY = 4'
                !Determining number of blocks and masks
                !CALL AllocIndex(block_size_qssa,mmxp,mmyp,mmzp,mia,miz,mja,mjz,mibcon, &
                !	     ngrids,dtlongn,N_DYN_CHEM)
	   
                CALL chem_qssa(mzp,mxp,myp,itchim,nr,mynum,ngrid,dtlt,basic_g(ngrid)%pp, &
           	 	       basic_g(ngrid)%pi0, basic_g(ngrid)%theta,basic_g(ngrid)%rv, &
           		       maxgrds,radiate_g(ngrid)%cosz,fast_JX_g(ngrid), &
           		       nspecies,nr_photo,weight,PhotoJMethod,chem1_g,index_g,nob, &
           		       maxblock_size)

                !Deallocating blocks data
                !CALL AllocIndex(block_size_qssa,mmxp,mmyp,mmzp,mia,miz,mja,mjz,mibcon, &
                !	     ngrids,dtlongn,N_DYN_CHEM)
	   
	   
             ELSEIF(CHEMISTRY == 2) THEN
	   	   
                CALL chem_ros(mzp,mxp,myp,maxgrds,ngrid,dtlt,basic_g(ngrid)%pp,basic_g(ngrid)%pi0, &
           	 	      basic_g(ngrid)%theta,basic_g(ngrid)%rv,basic_g(ngrid)%dn0, &
           		      radiate_g(ngrid)%cosz,micro_g(ngrid)%rcp,nspecies,nr,nr_photo,weight, &
           		      PhotojMethod,fast_JX_g(ngrid),maxnspecies,nspecies_chem_transported, &
           		      nspecies_chem_no_transported,transp_chem_index,no_transp_chem_index,chem1_g, &
           		      index_g,nob,maxblock_size,spack,spack_2d,get_non_zeros,&
			      N_DYN_CHEM,split_method)
             
             ELSEIF(CHEMISTRY == 3) THEN

                CALL chem_ros_dyndt(mzp,mxp,myp,maxgrds,ngrid,dtlt,basic_g(ngrid)%pp,basic_g(ngrid)%pi0, &
                  		    basic_g(ngrid)%theta,basic_g(ngrid)%rv,basic_g(ngrid)%dn0, &
                     	 	    radiate_g(ngrid)%cosz,micro_g(ngrid)%rcp,nspecies,nr,nr_photo,weight, &
                     		    PhotojMethod,maxnspecies,nspecies_chem_transported, &
                     		    nspecies_chem_no_transported,transp_chem_index,no_transp_chem_index, &
                     		    chem1_g,index_g,nob,maxblock_size,spack,spack_2d,atol,rtol,get_non_zeros,&
			            N_DYN_CHEM,split_method)
             ELSEIF(CHEMISTRY == 4) THEN              
  !print*,'84'; call flush(6)	  

                CALL chem_rodas3_dyndt(mzp,mxp,myp,maxgrds,ngrid,dtlt,basic_g(ngrid)%pp,basic_g(ngrid)%pi0, &
                     		       basic_g(ngrid)%theta,basic_g(ngrid)%rv,basic_g(ngrid)%dn0, &
                     		       radiate_g(ngrid)%cosz,micro_g(ngrid)%rcp,nspecies,nr,nr_photo,weight, &
                     		       PhotojMethod,maxnspecies,nspecies_chem_transported, &
                     		       nspecies_chem_no_transported,transp_chem_index,no_transp_chem_index, &
                     		       chem1_g,index_g,nob,maxblock_size,spack,spack_2d,atol,rtol,get_non_zeros,&
				       N_DYN_CHEM,split_method)
  !print*,'85'; call flush(6)	  
	     
             ENDIF

             !- set dyn tend to zero, only for parallel splliting method
             IF(split_method == 'PARALLEL' .AND. N_DYN_CHEM > 1) &
                  CALL chem_accum(ngrid,mzp,mxp,myp,ia,iz,ja,jz,2,nspecies_chem_transported,chem1_g, &
                                  transp_chem_index,N_DYN_CHEM)

          ENDIF

       ENDIF
       ! change MP 11/02/08
    ELSEIF(task == 4  .AND. CHEMISTRY  >= 1 .AND. CHEMISTRY_AQ >= 1 ) THEN
       CALL trans_gaq(mzp,mxp,myp,ia,iz,ja,jz)
     
    ELSEIF(task == 5  .AND. CHEMISTRY  >= 1 .AND. CHEMISTRY_AQ >= 1 ) THEN
       CALL trans_liq(mzp,mxp,myp,ia,iz,ja,jz,chem1_g,chem1aq_g,ngrid, &
                      nspeciesaq,ind_gas,chemic_g(ngrid)%coll,chemic_g(ngrid)%sedimr, &
                      micro_g(ngrid)%rcp,micro_g(ngrid)%rrp)
       ! change MP 11/02/08 - end
                
    ENDIF

  END SUBROUTINE chemistry_driver

!------------------------------------------------------------------------------------------

  SUBROUTINE include_tracer_lifetime(m1,m2,m3,spc_alloc_chem,nspecies,&
                                     transport,on,co,sc_p,sc_t)


    INTEGER , INTENT(IN)    :: m1
    INTEGER , INTENT(IN)    :: m2
    INTEGER , INTENT(IN)    :: m3
    INTEGER , INTENT(IN)    :: spc_alloc_chem(6,nspecies)
    INTEGER , INTENT(IN)    :: nspecies
    INTEGER , INTENT(IN)    :: transport
    INTEGER , INTENT(IN)    :: on
    INTEGER , INTENT(IN)    :: co
    REAL    , INTENT(IN)    :: sc_p(m1,m2,m3)
    REAL    , INTENT(INOUT) :: sc_t(m1*m2*m3)

    ! local declarations
    REAL,PARAMETER :: vm_CO_i   = 1./(30.*24.*3600.) ! 1/30 days
  
    ! - only for CO
    IF(spc_alloc_chem(transport,CO) == ON) THEN
   
       CALL accum_lifetime(m1*m2*m3     &
                           ,sc_t(1)     & ! tend
                           ,sc_p(1,1,1) & ! mix ratio
                           ,vm_CO_i   )   ! inv of lifetime (1/seconds)
    ENDIF
   
  END SUBROUTINE include_tracer_lifetime

!------------------------------------------------------------------------------------------

  SUBROUTINE accum_lifetime(m1,sc_t,sc_p,lf_inv)
       
    INTEGER , INTENT(IN)    :: m1
    REAL    , INTENT(INOUT) :: sc_t(m1)
    REAL    , INTENT(IN)    :: sc_p(m1)
    REAL    , INTENT(IN)    :: lf_inv

    INTEGER i

    DO i=1,m1
       sc_t(i) = sc_t(i) - lf_inv*sc_p(i)
    ENDDO
  END SUBROUTINE accum_lifetime

!------------------------------------------------------------------------------------------

  SUBROUTINE chem_accum(ngrid,m1,m2,m3,ia,iz,ja,jz,ntask,nspecies_chem_transported,chem1_g, &
                        transp_chem_index,N_DYN_CHEM)

     INTEGER         , INTENT(IN) :: ngrid
     INTEGER         , INTENT(IN) :: ntask
     INTEGER         , INTENT(IN) :: m1
     INTEGER         , INTENT(IN) :: m2
     INTEGER         , INTENT(IN) :: m3
     INTEGER         , INTENT(IN) :: ia
     INTEGER         , INTENT(IN) :: iz
     INTEGER         , INTENT(IN) :: ja
     INTEGER         , INTENT(IN) :: jz

    ! mem_chem1
    INTEGER          , INTENT(IN)    :: nspecies_chem_transported
    TYPE (chem1_vars), INTENT(INOUT) :: chem1_g(nspecies_chem_transported,ngrid)
    INTEGER          , INTENT(INOUT) :: transp_chem_index(nspecies_chem_transported)
    INTEGER          , INTENT(IN)    :: N_DYN_CHEM

   
    !- local var
    INTEGER ispc,n,ixyz,ntps,i,j,k
    REAL N_DYN_CHEM_i
    
    ntps = m1 * m2 * m3
       
    IF(ntask == 1) THEN
      
       N_DYN_CHEM_i= REAL(1./N_DYN_CHEM)
 	 
       DO ispc=1,nspecies_chem_transported 
    
          !- map the species to transported ones
          n=transp_chem_index(ispc)
      
          !- calculate the mean dynamic tendency for the entire chemistry timestep
 	  DO ixyz= 1,ntps
            chem1_g(n,ngrid)%sc_t_dyn(ixyz) = chem1_g(n,ngrid)%sc_t_dyn(ixyz) + &
                                              N_DYN_CHEM_i * chem1_g(n,ngrid)%sc_t(ixyz)
            !call accum(ntps, chem1_g(n,ngrid)%sc_t_dyn(1), chem1_g(n,ngrid)%sc_t(1) )
         ENDDO
 	  	  
      ENDDO
       
      !	do j=ja,jz ; do i=ia,iz ;do k=2,m1
      !	
      !
      !	!- Air temperature (K)  						  
      !	extra3d(9,ngrid)%d3(k,i,j)=extra3d(9,ngrid)%d3(k,i,j)+N_DYN_CHEM_i * &
      !				   basic_g(ngrid)%theta(k,i,j)  	* &
      !				( basic_g(ngrid)%pp (k,i,j) + basic_g(ngrid)%pi0(k,i,j) )/cp
      !
      !	!- transform from Exner function to pressure
      !	extra3d(10,ngrid)%d3(k,i,j)=extra3d(10,ngrid)%d3(k,i,j)+N_DYN_CHEM_i * &
      !				   p00* ( ( basic_g(ngrid)%pp (k,i,j) + basic_g(ngrid)%pi0(k,i,j) )/cp )**cpor
      !
      !	!- Water vapor 
      !	extra3d(11,ngrid)%d3(k,i,j)=extra3d(11,ngrid)%d3(k,i,j)+N_DYN_CHEM_i * &
      !				    basic_g(ngrid)% rv(k,i,j)
      !	enddo;enddo;enddo
       
   ELSE
      DO ispc=1,nspecies_chem_transported 
 	 
         n=transp_chem_index(ispc) !- map the species to transported ones
 	  
         !- set to zero arrays for  accumulation over the next chem timestep
         chem1_g(n,ngrid)%sc_t_dyn(1:ntps) = 0.
       
      ENDDO
      !	do j=ja,jz ; do i=ia,iz ;do k=2,m1
      !	extra3d(9 ,ngrid)%d3(k,i,j)=0.
      !	extra3d(10,ngrid)%d3(k,i,j)=0.
      !	extra3d(11,ngrid)%d3(k,i,j)=0.
      !	enddo;enddo;enddo
       
   ENDIF
 END  SUBROUTINE chem_accum
 !------------------------------------------------------------------
  
  
END MODULE module_chemistry_driver
  
  
  

































SUBROUTINE latset_tracer(m1,m2,m3,ia,iz,ja,jz,ibcon,ap,uc,vc,dxu, &
     dxm,dyv,dym)
  USE mem_grid, ONLY: dtlt,ibnd,jdim,jbnd,nz
  USE mem_scratch, ONLY: vctr17,vctr18
  IMPLICIT NONE

  INTEGER :: m1,m2,m3,ia,iz,ja,jz,ibcon,i,j,k,lbw,lbe,lbs,lbn
  REAL :: thresh,dtlx,c1,dxr,dyr
  REAL, DIMENSION(m1,m2,m3) :: ap,uc,vc
  REAL, DIMENSION(m2,m3) :: dxu,dxm,dyv,dym
  !character(len=*) :: vnam

  IF (IAND(ibcon,1) .GT. 0) lbw = ia - 1
  IF (IAND(ibcon,2) .GT. 0) lbe = iz + 1
  IF (IAND(ibcon,4) .GT. 0) lbs = ja - 1
  IF (IAND(ibcon,8) .GT. 0) lbn = jz + 1

  thresh = 0.
  dtlx = dtlt


  IF (ibnd .NE. 4) THEN

     ! Western boundary for lsflg = 2  == constant inflow, radiative b.c. outflow
     ! Veja que o campo ap(k,i,j) na borda so e' modificado se houver outflow, isto e'
     ! em regioes de inflow o valor permanece constante e igual ao valor inicial (t=0).
     ! Para introduzir uma especie de "nudging" na condicao de contorno, modifique o
     ! vetor ap(k,lbw,j) para a 'western boundary', por exemplo. O mesmo para as outras
     ! bordas.

     IF (IAND(ibcon,1) .GT. 0) THEN
	DO j = 1,m3
	   dxr = dxu(ia,j) / dxu(lbw,j)
	   c1 = dtlx * dxu(lbw,j)
	   DO k = 1,m1
	      vctr17(k) = -c1 * uc(k,lbw,j)
	      vctr18(k) = ap(k,ia,j) + dxr * (ap(k,ia,j) - ap(k,ia+1,j))
	   ENDDO
	   DO k = 1,m1
	      IF (vctr17(k) .GE. thresh) ap(k,lbw,j) = vctr18(k)
	   ENDDO
	ENDDO
     ENDIF

     !     Eastern Boundary for LSFLG =  2

     IF (IAND(ibcon,2) .GT. 0) THEN
	DO j = 1,m3
	   dxr = dxu(iz-1,j) / dxu(iz,j)
	   c1 = dtlx * dxu(iz,j)
	   DO k = 1,m1
	      vctr17(k) = c1 * uc(k,iz,j)
	      vctr18(k) = ap(k,iz,j) + dxr * (ap(k,iz,j) - ap(k,iz-1,j))
	   ENDDO
	   DO k = 1,m1
	      IF (vctr17(k) .GE. thresh)  ap(k,lbe,j) = vctr18(k)
	   ENDDO
	ENDDO
     ENDIF
  ENDIF

  IF(jdim.EQ.1.AND.jbnd.NE.4)THEN

     !     Southern boundary for LSFLG  2

     IF (IAND(ibcon,4) .GT. 0) THEN
	DO i = 1,m2
	   dyr = dyv(i,ja) / dyv(i,lbs)
	   c1 = dtlx * dyv(i,lbs)
!srf - fix from rams 60
!             do k = 1,nz
              DO k = 1,m1
	      vctr17(k) = -c1 * vc(k,i,lbs)
	      vctr18(k) = ap(k,i,ja) + dyr * (ap(k,i,ja) - ap(k,i,ja+1))
	   ENDDO
	   DO k = 1,m1
	      IF (vctr17(k) .GE. thresh) ap(k,i,lbs) = vctr18(k)
	   ENDDO
	ENDDO
     ENDIF

     !     Northern Boundary for LSFLG =  2

     IF (IAND(ibcon,8) .GT. 0) THEN
	DO i = 1,m2
	   dyr = dyv(i,jz-1) / dyv(i,jz)
	   c1 = dtlx * dyv(i,jz)
	   DO k = 1,m1
	      vctr17(k) = c1 * vc(k,i,jz)
	      vctr18(k) = ap(k,i,jz) + dyr * (ap(k,i,jz) - ap(k,i,jz-1))
	   ENDDO
	   DO k = 1,m1
	      IF (vctr17(k) .GE. thresh) ap(k,i,lbn) = vctr18(k)
	   ENDDO
	ENDDO
     ENDIF
  ENDIF

END SUBROUTINE latset_tracer
!-----------------------------------------------------------------------

!------------------------------------------------------------------
  SUBROUTINE initial_condition(ng,m1,m2,m3) 
    USE chem1_list, ONLY : NX=>nspecies,init_ajust,CO 
    USE mem_chem1, ONLY: chem1_g,chem1_src_g,CHEM_ASSIM
    USE mem_basic         , ONLY: basic_g
    USE aer1_list
    USE mem_aer1
    IMPLICIT NONE
    INTEGER :: ng,m1,m2,m3,i,j,k,ispc
    REAL,PARAMETER :: fcu =1.e+9 !=> ppbm 
   
   
    IF(CHEM_ASSIM > 0 ) RETURN
   
    !-attention: values must be in ppbv (check if this correct, should 
    !                                    not be in ppbm? SRF)
    !!print*,'============================= homogeneous initial. ====='
    !!print*,'CO= ',chem1_g(CO ,ng)%sc_p(1:m1,72,68)*fcu
    !!print*,'========================================================'
    DO ispc=1,NX
!- for homogeneous 100 pppv everywhere
!    chem1_g(ispc,ng)%sc_p(:,:,:)=1.*100.*(28/28.96)
!- for a box of cte mixing ratio 
!     DO j=30,50
!      DO i=30,50
!       DO k=10,12
!         chem1_g(ispc,ng)%sc_p(k,i,j)=100.*(28/28.96)
!    ENDDO;ENDDO;ENDDO;ENDDO

    DO j=1,m3
      DO i=1,m2
       DO k=1,m1
!         chem1_g(ispc,ng)%sc_p(k,i,j)=chem1_g(ispc ,ng)%sc_p(k,72,68)
         chem1_g(ispc,ng)%sc_p(k,i,j)=float(j)
    ENDDO;ENDDO;ENDDO;ENDDO
    
!- to apply init_ajust factor
    DO ispc=1,nx
      chem1_g(ispc,ng)%sc_p(:,:,:)=init_ajust(ispc)*chem1_g(ispc ,ng)%sc_p(:,:,:)
    ENDDO
  
  !stop 33
    RETURN
    
    !
    !do i=1,m2
    ! do j=1,m3
    ! chem1_g(O3 ,ng)%sc_p(1:13,i,j)=10. + i*j/100
    !  chem1_g(no ,ng)%sc_p(1:20,i,j)=20. + i*j**2/1000
    !  chem1_g(no2 ,ng)%sc_p(1:10,i,j)=1. + (i**2) * j*2/1000
    !enddo
    ! enddo
    !chem1_g(H2O,ng)%sc_p(:,:,:)=basic_g(ng)%RV(:,:,:)*fcu
    
    RETURN
    
   ! chem1_g(CO,ng)%sc_p(:,:,:) = 60.
    !chem1_g(O3,ng)%sc_p(:,:,:) = 50.
    !chem1_g(O3,ng)%sc_p(25:,:,:) = 100.
    
    !chem1_g(H2O2,ng)%sc_p(1:13,:,:) = 9.
    
    aer1_g(accum,bburn,ng)%sc_p(1:20,3:m2-3,3:m3-3) = 10.
    aer1_g(coarse,bburn,ng)%sc_p(1:10,3:m2-3,3:m3-3) = 10.
    aer1_g(nucle,urban,ng)%sc_p(1:10,3:m2-3,3:m3-3) = 10.
    aer1_g(accum,urban,ng)%sc_p(1:10,3:m2-3,3:m3-3) = 10.
    
    aer1_g(accum,bburn,ng)%sc_p(17:19,3:m2-3,3:m3-3) = 8.
    aer1_g(coarse,bburn,ng)%sc_p(17:19,3:m2-3,3:m3-3) = 8.
    aer1_g(nucle,urban,ng)%sc_p(17:19,3:m2-3,3:m3-3) = 8.
    aer1_g(accum,urban,ng)%sc_p(17:19,3:m2-3,3:m3-3) = 8.


  END SUBROUTINE initial_condition


  	
